Untitled
Introduction
Host plant genotype and phenotype influence arthropod community composition. They influence fitness, can mediate competition, predation and host-parasite interactions. Plant size, leaf chemistry, defense compounds, trichome density etc all affect species in different ways. However, it is hard to pin down general principles in how host plant characteristics affect herbivores across species Insect traits can be used to provide a mechanistic connection between species and environment. Insect traits have increasingly been used to understand insect response to characteristics of landscape and understand community composition.
I set out with the following objectives:
(i) How do herbivore communities respond to phenotypic variation?
(ii) Which plant traits account for herbivore community responses to host-plant phenotype?
(iii) Which insect traits account for herbivore community responses to host-plant phenotype
Results
Some context: Richness across hybrid sites
(i) How do herbivore communities respond to phenotypic variation?
p<-metaMDS(hem.both, distance="bray", k=2, trymax=100)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2614694
## Run 1 stress 0.2591777
## ... New best solution
## ... Procrustes: rmse 0.07578497 max resid 0.2182768
## Run 2 stress 0.2643395
## Run 3 stress 0.2711564
## Run 4 stress 0.2591012
## ... New best solution
## ... Procrustes: rmse 0.004049879 max resid 0.02030985
## Run 5 stress 0.2680823
## Run 6 stress 0.2669241
## Run 7 stress 0.2676229
## Run 8 stress 0.2668678
## Run 9 stress 0.2686223
## Run 10 stress 0.264573
## Run 11 stress 0.2650065
## Run 12 stress 0.2764339
## Run 13 stress 0.3046596
## Run 14 stress 0.2603386
## Run 15 stress 0.2646121
## Run 16 stress 0.2642043
## Run 17 stress 0.260578
## Run 18 stress 0.264796
## Run 19 stress 0.2646198
## Run 20 stress 0.2642683
## Run 21 stress 0.2708115
## Run 22 stress 0.261162
## Run 23 stress 0.2654999
## Run 24 stress 0.2630344
## Run 25 stress 0.2657158
## Run 26 stress 0.2587648
## ... New best solution
## ... Procrustes: rmse 0.01136149 max resid 0.04027083
## Run 27 stress 0.2617533
## Run 28 stress 0.263632
## Run 29 stress 0.2749018
## Run 30 stress 0.2654376
## Run 31 stress 0.2621657
## Run 32 stress 0.2667588
## Run 33 stress 0.265062
## Run 34 stress 0.2766734
## Run 35 stress 0.2720742
## Run 36 stress 0.2680699
## Run 37 stress 0.2671735
## Run 38 stress 0.2703416
## Run 39 stress 0.2651467
## Run 40 stress 0.2658195
## Run 41 stress 0.2629638
## Run 42 stress 0.2648128
## Run 43 stress 0.2630296
## Run 44 stress 0.2613272
## Run 45 stress 0.2648603
## Run 46 stress 0.2618145
## Run 47 stress 0.27441
## Run 48 stress 0.267981
## Run 49 stress 0.266718
## Run 50 stress 0.2606572
## Run 51 stress 0.2696575
## Run 52 stress 0.2667685
## Run 53 stress 0.2636413
## Run 54 stress 0.2686025
## Run 55 stress 0.2697989
## Run 56 stress 0.2589156
## ... Procrustes: rmse 0.009479284 max resid 0.04055719
## Run 57 stress 0.2648442
## Run 58 stress 0.2603029
## Run 59 stress 0.2696071
## Run 60 stress 0.2615379
## Run 61 stress 0.2652593
## Run 62 stress 0.2659076
## Run 63 stress 0.2666821
## Run 64 stress 0.2682568
## Run 65 stress 0.2668895
## Run 66 stress 0.2711279
## Run 67 stress 0.2657382
## Run 68 stress 0.2654936
## Run 69 stress 0.2676306
## Run 70 stress 0.2818346
## Run 71 stress 0.2822666
## Run 72 stress 0.2655272
## Run 73 stress 0.2660869
## Run 74 stress 0.264887
## Run 75 stress 0.2619298
## Run 76 stress 0.2634607
## Run 77 stress 0.2650585
## Run 78 stress 0.2708242
## Run 79 stress 0.2644942
## Run 80 stress 0.2731182
## Run 81 stress 0.2611793
## Run 82 stress 0.2629705
## Run 83 stress 0.2653446
## Run 84 stress 0.2772961
## Run 85 stress 0.2698972
## Run 86 stress 0.2651032
## Run 87 stress 0.2745211
## Run 88 stress 0.2655207
## Run 89 stress 0.2652121
## Run 90 stress 0.263763
## Run 91 stress 0.2619298
## Run 92 stress 0.2654345
## Run 93 stress 0.2731397
## Run 94 stress 0.2630328
## Run 95 stress 0.2627176
## Run 96 stress 0.2668713
## Run 97 stress 0.2614056
## Run 98 stress 0.2603036
## Run 99 stress 0.265145
## Run 100 stress 0.2633985
## *** No convergence -- monoMDS stopping criteria:
## 100: stress ratio > sratmaxordiplot(p,type="n")
orditorp(p,display="species",col="red",air=0.01)
orditorp(p,display="sites",cex=1.25,air=0.01)
ordiplot(p,type="n")
orditorp(p,display="species",col="red",air=0.01)
orditorp(p,display="sites",pch=19, col="blue", label=F)
ordiellipse(p,groups=env.both$morphotype,draw="polygon",col="grey90",label=T, kind="se")
Are communities distinctly different between different morphotype trees within the same site?
hem.ero<-hem[env$site=="ERO",]
hem.ali<-hem[env$site=="ALI",]
hem.kaiy<-hem[env$site=="KAIY",]
hem.olaa<-hem[env$site=="OLAA",]
hem.kaio<-hem[env$site=="KAIO",]
#par(mfrow=c(2,2))
moment<-metaMDS(hem.ero, k=3)## Wisconsin double standardization
## Run 0 stress 0.06827439
## Run 1 stress 0.06827615
## ... Procrustes: rmse 0.0004506555 max resid 0.0007439575
## ... Similar to previous best
## Run 2 stress 0.06827984
## ... Procrustes: rmse 0.0008459453 max resid 0.001391683
## ... Similar to previous best
## Run 3 stress 0.068273
## ... New best solution
## ... Procrustes: rmse 0.0003991526 max resid 0.0007228103
## ... Similar to previous best
## Run 4 stress 0.08827491
## Run 5 stress 0.0822194
## Run 6 stress 0.08132868
## Run 7 stress 0.08854191
## Run 8 stress 0.06827149
## ... New best solution
## ... Procrustes: rmse 0.003514993 max resid 0.006012255
## ... Similar to previous best
## Run 9 stress 0.1413167
## Run 10 stress 0.0882751
## Run 11 stress 0.06827118
## ... New best solution
## ... Procrustes: rmse 0.002863892 max resid 0.005034434
## ... Similar to previous best
## Run 12 stress 0.08132942
## Run 13 stress 0.06827011
## ... New best solution
## ... Procrustes: rmse 0.0005564563 max resid 0.001017243
## ... Similar to previous best
## Run 14 stress 0.08133052
## Run 15 stress 0.08827486
## Run 16 stress 0.2421808
## Run 17 stress 0.06828125
## ... Procrustes: rmse 0.00228561 max resid 0.003866028
## ... Similar to previous best
## Run 18 stress 0.06827822
## ... Procrustes: rmse 0.002354057 max resid 0.003952792
## ... Similar to previous best
## Run 19 stress 0.06827094
## ... Procrustes: rmse 0.002173734 max resid 0.003670135
## ... Similar to previous best
## Run 20 stress 0.06827642
## ... Procrustes: rmse 0.001993287 max resid 0.003361344
## ... Similar to previous best
## *** Solution reachedmds.fig <- ordiplot(moment, type = "none", main="Escape road old flow")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="ERO",]$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env[env$site=="ERO",]$morphotype == "P")
Lab<-c( "Hybrid", "Pubescent")
color<-c("blue", "red")
legend( "topright", legend=Lab,pch=c(15,19),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)
moment<-metaMDS(hem.kaiy, k=3)## Wisconsin double standardization
## Run 0 stress 0.06676995
## Run 1 stress 0.06676756
## ... New best solution
## ... Procrustes: rmse 0.0004757995 max resid 0.000815727
## ... Similar to previous best
## Run 2 stress 0.08024762
## Run 3 stress 0.0802441
## Run 4 stress 0.08845382
## Run 5 stress 0.08024536
## Run 6 stress 0.06676711
## ... New best solution
## ... Procrustes: rmse 0.004773203 max resid 0.008547457
## ... Similar to previous best
## Run 7 stress 0.08627771
## Run 8 stress 0.08024788
## Run 9 stress 0.08627932
## Run 10 stress 0.0667668
## ... New best solution
## ... Procrustes: rmse 0.0001037532 max resid 0.0002291509
## ... Similar to previous best
## Run 11 stress 0.08848844
## Run 12 stress 0.08847556
## Run 13 stress 0.066772
## ... Procrustes: rmse 0.0008988719 max resid 0.001495977
## ... Similar to previous best
## Run 14 stress 0.08627253
## Run 15 stress 0.06676287
## ... New best solution
## ... Procrustes: rmse 0.001662798 max resid 0.002918259
## ... Similar to previous best
## Run 16 stress 0.06676659
## ... Procrustes: rmse 0.001537349 max resid 0.002704864
## ... Similar to previous best
## Run 17 stress 0.06677093
## ... Procrustes: rmse 0.003762415 max resid 0.006646404
## ... Similar to previous best
## Run 18 stress 0.06676934
## ... Procrustes: rmse 0.003467202 max resid 0.006118675
## ... Similar to previous best
## Run 19 stress 0.06676371
## ... Procrustes: rmse 0.0007865925 max resid 0.001268778
## ... Similar to previous best
## Run 20 stress 0.06676702
## ... Procrustes: rmse 0.002951224 max resid 0.005323313
## ... Similar to previous best
## *** Solution reachedmds.fig <- ordiplot(moment, type = "none", main="Kaiholena Young flow")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="KAIY",]$morphotype == "G")
points(mds.fig, "sites", pch = 19, col = "red", select = env[env$site=="KAIY",]$morphotype == "P")
Lab<-c("Glabrous", "Pubescent")
color<-c("blue", "red")
legend( "topright", legend=Lab,pch=c(22,19),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)
moment<-metaMDS(hem.ali, k=3)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.0713509
## Run 1 stress 0.07135184
## ... Procrustes: rmse 0.0003176795 max resid 0.0006529927
## ... Similar to previous best
## Run 2 stress 0.07135107
## ... Procrustes: rmse 0.0003125567 max resid 0.0006357934
## ... Similar to previous best
## Run 3 stress 0.09581055
## Run 4 stress 0.09581243
## Run 5 stress 0.07135091
## ... Procrustes: rmse 0.0002325267 max resid 0.0004489137
## ... Similar to previous best
## Run 6 stress 0.09581176
## Run 7 stress 0.09581325
## Run 8 stress 0.07135229
## ... Procrustes: rmse 0.0005044818 max resid 0.001098918
## ... Similar to previous best
## Run 9 stress 0.07135072
## ... New best solution
## ... Procrustes: rmse 0.0002041099 max resid 0.0004146758
## ... Similar to previous best
## Run 10 stress 0.07135078
## ... Procrustes: rmse 0.0001761285 max resid 0.0003919795
## ... Similar to previous best
## Run 11 stress 0.07135076
## ... Procrustes: rmse 6.081179e-05 max resid 0.000105218
## ... Similar to previous best
## Run 12 stress 0.07135164
## ... Procrustes: rmse 0.0003903769 max resid 0.0009003304
## ... Similar to previous best
## Run 13 stress 0.07135088
## ... Procrustes: rmse 0.0001376946 max resid 0.0002526987
## ... Similar to previous best
## Run 14 stress 0.07135159
## ... Procrustes: rmse 0.000382806 max resid 0.0008377661
## ... Similar to previous best
## Run 15 stress 0.07135072
## ... Procrustes: rmse 7.440481e-05 max resid 0.0001546999
## ... Similar to previous best
## Run 16 stress 0.07135074
## ... Procrustes: rmse 0.00012075 max resid 0.0002683365
## ... Similar to previous best
## Run 17 stress 0.0713508
## ... Procrustes: rmse 5.576722e-05 max resid 9.502438e-05
## ... Similar to previous best
## Run 18 stress 0.07135075
## ... Procrustes: rmse 4.769923e-05 max resid 8.960165e-05
## ... Similar to previous best
## Run 19 stress 0.07135082
## ... Procrustes: rmse 0.0001276776 max resid 0.0002492422
## ... Similar to previous best
## Run 20 stress 0.07135201
## ... Procrustes: rmse 0.0004838851 max resid 0.001138584
## ... Similar to previous best
## *** Solution reachedmds.fig <- ordiplot(moment, type = "none", main="Alili Springs")
text(mds.fig, "species",cex=0.5)
points(mds.fig, "sites", pch = 19, col = "green", select = env[env$site=="ALI",]$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="ALI",]$morphotype == "H")
Lab<-c("G", "H")
color<-c("green", "blue")
legend( "topright", legend=Lab,pch=c(19,15),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)
moment<-metaMDS(hem.olaa, k=3)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.06142282
## Run 1 stress 0.06142317
## ... Procrustes: rmse 0.0001459225 max resid 0.0002542925
## ... Similar to previous best
## Run 2 stress 0.06142369
## ... Procrustes: rmse 0.001569287 max resid 0.002309513
## ... Similar to previous best
## Run 3 stress 0.06142404
## ... Procrustes: rmse 0.001889031 max resid 0.003056083
## ... Similar to previous best
## Run 4 stress 0.06142421
## ... Procrustes: rmse 0.0004712779 max resid 0.0007472414
## ... Similar to previous best
## Run 5 stress 0.06142531
## ... Procrustes: rmse 0.0006532095 max resid 0.0009522755
## ... Similar to previous best
## Run 6 stress 0.06142253
## ... New best solution
## ... Procrustes: rmse 0.0004460018 max resid 0.0006143714
## ... Similar to previous best
## Run 7 stress 0.06177549
## ... Procrustes: rmse 0.01106618 max resid 0.01887957
## Run 8 stress 0.06142506
## ... Procrustes: rmse 0.001144175 max resid 0.001729107
## ... Similar to previous best
## Run 9 stress 0.06142314
## ... Procrustes: rmse 0.0004750837 max resid 0.0007324986
## ... Similar to previous best
## Run 10 stress 0.06142336
## ... Procrustes: rmse 0.0006656519 max resid 0.0009760464
## ... Similar to previous best
## Run 11 stress 0.06142309
## ... Procrustes: rmse 0.0006042221 max resid 0.0008453728
## ... Similar to previous best
## Run 12 stress 0.06652333
## Run 13 stress 0.06142247
## ... New best solution
## ... Procrustes: rmse 0.0003501391 max resid 0.0006710117
## ... Similar to previous best
## Run 14 stress 0.06652333
## Run 15 stress 0.06142431
## ... Procrustes: rmse 0.0008082026 max resid 0.001309491
## ... Similar to previous best
## Run 16 stress 0.06142235
## ... New best solution
## ... Procrustes: rmse 0.0001497941 max resid 0.0002665819
## ... Similar to previous best
## Run 17 stress 0.06652356
## Run 18 stress 0.06652488
## Run 19 stress 0.06142374
## ... Procrustes: rmse 0.0005622079 max resid 0.000919252
## ... Similar to previous best
## Run 20 stress 0.06142262
## ... Procrustes: rmse 0.0001982676 max resid 0.0002924331
## ... Similar to previous best
## *** Solution reachedmds.fig <- ordiplot(moment, type = "none", main="HAVO Olaa")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 19, col = "green", select = env[env$site=="OLAA",]$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="OLAA",]$morphotype == "H")
Lab<-c("Glabrous", "Hybrid")
color<-c("green", "blue" )
legend( "topright", legend=Lab,pch=c(19,15),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)
# experiment with adding data kaiholena old- only 6 data points
moment<-metaMDS(hem.kaio, k=3)## Wisconsin double standardization
## Run 0 stress 0
## Run 1 stress 0
## ... Procrustes: rmse 0.1099498 max resid 0.1686382
## Run 2 stress 9.310127e-05
## ... Procrustes: rmse 0.08995409 max resid 0.1557758
## Run 3 stress 8.792988e-05
## ... Procrustes: rmse 0.1389952 max resid 0.1699884
## Run 4 stress 0
## ... Procrustes: rmse 0.0510628 max resid 0.07329546
## Run 5 stress 0
## ... Procrustes: rmse 0.1677124 max resid 0.2726518
## Run 6 stress 0
## ... Procrustes: rmse 0.1996001 max resid 0.2891551
## Run 7 stress 6.345724e-05
## ... Procrustes: rmse 0.2144208 max resid 0.2853603
## Run 8 stress 0
## ... Procrustes: rmse 0.1989755 max resid 0.3117667
## Run 9 stress 9.697878e-05
## ... Procrustes: rmse 0.1731629 max resid 0.2866622
## Run 10 stress 0
## ... Procrustes: rmse 0.1202304 max resid 0.1895108
## Run 11 stress 8.569755e-05
## ... Procrustes: rmse 0.1706865 max resid 0.2845998
## Run 12 stress 9.006129e-05
## ... Procrustes: rmse 0.1727758 max resid 0.2864063
## Run 13 stress 8.904747e-05
## ... Procrustes: rmse 0.2276194 max resid 0.3055194
## Run 14 stress 0
## ... Procrustes: rmse 0.1480634 max resid 0.2556493
## Run 15 stress 9.791498e-05
## ... Procrustes: rmse 0.2186505 max resid 0.3068564
## Run 16 stress 8.368217e-05
## ... Procrustes: rmse 0.1627448 max resid 0.2784951
## Run 17 stress 0
## ... Procrustes: rmse 0.1701545 max resid 0.2513082
## Run 18 stress 9.282268e-05
## ... Procrustes: rmse 0.21155 max resid 0.3065631
## Run 19 stress 0
## ... Procrustes: rmse 0.1465121 max resid 0.2101937
## Run 20 stress 0
## ... Procrustes: rmse 0.1146521 max resid 0.1938358
## *** No convergence -- monoMDS stopping criteria:
## 20: stress < smin## Warning in metaMDS(hem.kaio, k = 3): stress is (nearly) zero: you may have
## insufficient data## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below thresholdmds.fig <- ordiplot(moment, type = "none", main="KAIO limited data")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env[env$site=="KAIO",]$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env[env$site=="KAIO",]$morphotype == "P")
points(mds.fig, "sites", pch = 19, col = "green", select = env[env$site=="KAIO",]$morphotype == "G")
Lab<-c( "Hybrid", "Pubescent", "Glabrous")
color<-c("blue", "red", "green")
legend( "topright", legend=Lab,pch=c(15,19),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)
# explore adding kaio to kaiy.
hem.kai.all<-rbind(hem.kaiy,hem.kaio)
env.kai.all<-rbind(env[env$site=="KAIY",],env[env$site=="KAIO",])
moment<-metaMDS(hem.kai.all, k=3)## Wisconsin double standardization
## Run 0 stress 0.1432029
## Run 1 stress 0.1400188
## ... New best solution
## ... Procrustes: rmse 0.1367607 max resid 0.328315
## Run 2 stress 0.1450392
## Run 3 stress 0.1400853
## ... Procrustes: rmse 0.1416568 max resid 0.4080269
## Run 4 stress 0.1349444
## ... New best solution
## ... Procrustes: rmse 0.1182893 max resid 0.2461261
## Run 5 stress 0.1349426
## ... New best solution
## ... Procrustes: rmse 0.0007411143 max resid 0.002142872
## ... Similar to previous best
## Run 6 stress 0.1369001
## Run 7 stress 0.1437608
## Run 8 stress 0.1465123
## Run 9 stress 0.1401
## Run 10 stress 0.1450611
## Run 11 stress 0.1386605
## Run 12 stress 0.1369041
## Run 13 stress 0.1349427
## ... Procrustes: rmse 0.00028015 max resid 0.0005560841
## ... Similar to previous best
## Run 14 stress 0.1400167
## Run 15 stress 0.1400186
## Run 16 stress 0.1465886
## Run 17 stress 0.1349427
## ... Procrustes: rmse 0.0006100189 max resid 0.001454125
## ... Similar to previous best
## Run 18 stress 0.1400879
## Run 19 stress 0.1400927
## Run 20 stress 0.1400951
## *** Solution reachedmds.fig <- ordiplot(moment, type = "none", main="KAI all")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env.kai.all$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env.kai.all$morphotype == "P")
points(mds.fig, "sites", pch = 19, col = "green", select = env.kai.all$morphotype == "G")
adonis(hem.kai.all~env.kai.all$morphotype, by="margin")##
## Call:
## adonis(formula = hem.kai.all ~ env.kai.all$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## env.kai.all$morphotype 2 1.4422 0.72111 2.5311 0.24034 0.001 ***
## Residuals 16 4.5584 0.28490 0.75966
## Total 18 6.0006 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(hem.ali~env[env$site=="ALI",]$morphotype, by="margin")##
## Call:
## adonis(formula = hem.ali ~ env[env$site == "ALI", ]$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2
## env[env$site == "ALI", ]$morphotype 1 0.35248 0.35248 1.4706 0.12821
## Residuals 10 2.39678 0.23968 0.87179
## Total 11 2.74926 1.00000
## Pr(>F)
## env[env$site == "ALI", ]$morphotype 0.174
## Residuals
## Totaladonis(hem.kaiy~env[env$site=="KAIY",]$morphotype, by="margin")##
## Call:
## adonis(formula = hem.kaiy ~ env[env$site == "KAIY", ]$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2
## env[env$site == "KAIY", ]$morphotype 1 1.0913 1.09131 4.5097 0.29077
## Residuals 11 2.6619 0.24199 0.70923
## Total 12 3.7532 1.00000
## Pr(>F)
## env[env$site == "KAIY", ]$morphotype 0.002 **
## Residuals
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(hem.ero~env[env$site=="ERO",]$morphotype, by="margin")##
## Call:
## adonis(formula = hem.ero ~ env[env$site == "ERO", ]$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2
## env[env$site == "ERO", ]$morphotype 1 0.4346 0.43459 1.3671 0.12027
## Residuals 10 3.1789 0.31789 0.87973
## Total 11 3.6135 1.00000
## Pr(>F)
## env[env$site == "ERO", ]$morphotype 0.158
## Residuals
## Totaladonis(hem.olaa~env[env$site=="OLAA",]$morphotype, by="margin")##
## Call:
## adonis(formula = hem.olaa ~ env[env$site == "OLAA", ]$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2
## env[env$site == "OLAA", ]$morphotype 1 0.7505 0.75048 2.5677 0.20431
## Residuals 10 2.9228 0.29228 0.79569
## Total 11 3.6733 1.00000
## Pr(>F)
## env[env$site == "OLAA", ]$morphotype 0.045 *
## Residuals
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There are differences in species richness as well as species composition across sites. Diversity is highest at the high productivity sites midway the Big Island chronosequence. Species composition varies between sites, though many species are shared in common.
(ii) Which plant traits account for herbivore community responses to host-plant phenotype?
# partition variance between environment and spatial
#test<-rda(decostand(trial, "hel"), leaftraits[,c(1:4,6)], xy.2014[,2:3])
#test<-varpart(vegdist(hem), env$z.age, xy[,2:3])
#capscale(vegdist(trial)~ leaftraits[,1:4]+xy[,2:3])
#test
#plot(test)
test<-varpart(vegdist(decostand(hem.2014, "hel")), leaftraits[,c(1,2,3,4,6)], xy.2014[,2:3])
test<-rda(decostand(hem.2014, "hel"), leaftraits[,c(1,2,3,4,6)], xy.2014[,2:3])
test.wout<-rda(decostand(hem.2014, "hel"), leaftraits[,c(1,2,3,4,6)])
test.spat<-rda(decostand(hem.2014, "hel"), xy.2014[,2:3])
anova(test)## Permutation test for rda under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model: rda(X = decostand(hem.2014, "hel"), Y = leaftraits[, c(1, 2, 3, 4, 6)], Z = xy.2014[, 2:3])
## Df Variance F Pr(>F)
## Model 6 0.22081 1.7192 0.008 **
## Residual 15 0.32109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(hem.2014~leaftraits$percN+leaftraits$meanWatercontent+leaftraits$meanSLA+leaftraits$percP+xy.2014$long+env.2014$morphotype+xy.2014$lat, by="margin")##
## Call:
## adonis(formula = hem.2014 ~ leaftraits$percN + leaftraits$meanWatercontent + leaftraits$meanSLA + leaftraits$percP + xy.2014$long + env.2014$morphotype + xy.2014$lat, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## leaftraits$percN 1 1.1225 1.12250 4.4231 0.13512 0.001
## leaftraits$meanWatercontent 1 0.4434 0.44340 1.7472 0.05338 0.054
## leaftraits$meanSLA 1 0.3591 0.35912 1.4151 0.04323 0.161
## leaftraits$percP 1 0.9557 0.95568 3.7658 0.11504 0.001
## xy.2014$long 1 0.5235 0.52347 2.0627 0.06301 0.025
## env.2014$morphotype 2 0.6836 0.34179 1.3468 0.08229 0.136
## xy.2014$lat 1 0.4127 0.41273 1.6263 0.04968 0.086
## Residuals 15 3.8067 0.25378 0.45824
## Total 23 8.3072 1.00000
##
## leaftraits$percN ***
## leaftraits$meanWatercontent .
## leaftraits$meanSLA
## leaftraits$percP ***
## xy.2014$long *
## env.2014$morphotype
## xy.2014$lat .
## Residuals
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(hem.2014~xy.2014$lat+xy.2014$long+leaftraits$percN+leaftraits$percP+leaftraits$meanSLA+leaftraits$meanWatercontent+env.2014$morphotype, by="margin")##
## Call:
## adonis(formula = hem.2014 ~ xy.2014$lat + xy.2014$long + leaftraits$percN + leaftraits$percP + leaftraits$meanSLA + leaftraits$meanWatercontent + env.2014$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## xy.2014$lat 1 1.1255 1.12547 4.4349 0.13548 0.001
## xy.2014$long 1 0.7394 0.73945 2.9137 0.08901 0.005
## leaftraits$percN 1 0.9300 0.93003 3.6647 0.11196 0.001
## leaftraits$percP 1 0.3141 0.31412 1.2378 0.03781 0.275
## leaftraits$meanSLA 1 0.3360 0.33602 1.3241 0.04045 0.195
## leaftraits$meanWatercontent 1 0.4945 0.49448 1.9485 0.05952 0.034
## env.2014$morphotype 2 0.5609 0.28045 1.1051 0.06752 0.341
## Residuals 15 3.8067 0.25378 0.45824
## Total 23 8.3072 1.00000
##
## xy.2014$lat ***
## xy.2014$long **
## leaftraits$percN ***
## leaftraits$percP
## leaftraits$meanSLA
## leaftraits$meanWatercontent *
## env.2014$morphotype
## Residuals
## Total
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# rda/varpart: when testing only 2014 hybrid zone sites from 2014, xy coordinates explain some. better explained still by leaftraits
# adonis seems to indicate latitude is important predictor- sig even after all other leaftrait variables
# while almost always sig (for 170 data points and 24), R2 is less than 0.1 in all cases
adonis(hem.2014~leaftraits$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$meanSLA+leaftraits$percP+scale(xy.2014$lat,scale=T,center=T), by="margin")##
## Call:
## adonis(formula = hem.2014 ~ leaftraits$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$meanSLA + leaftraits$percP + scale(xy.2014$lat, scale = T, center = T), by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model
## leaftraits$morphotype 2 1.3702 0.68512 2.6179
## leaftraits$percN 1 0.6270 0.62697 2.3957
## leaftraits$meanWatercontent 1 0.5596 0.55960 2.1383
## leaftraits$meanSLA 1 0.3781 0.37811 1.4448
## leaftraits$percP 1 0.6230 0.62301 2.3806
## scale(xy.2014$lat, scale = T, center = T) 1 0.5619 0.56193 2.1472
## Residuals 16 4.1873 0.26171
## Total 23 8.3072
## R2 Pr(>F)
## leaftraits$morphotype 0.16495 0.001 ***
## leaftraits$percN 0.07547 0.006 **
## leaftraits$meanWatercontent 0.06736 0.032 *
## leaftraits$meanSLA 0.04552 0.153
## leaftraits$percP 0.07500 0.012 *
## scale(xy.2014$lat, scale = T, center = T) 0.06764 0.036 *
## Residuals 0.50406
## Total 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1adonis(hem.2014~env.2014$tree.height+leaftraits$morphotype+leaftraits$meanSLA+leaftraits$percN+env.2014$z.age+xy.2014$lat+xy.2014$long+env.2014$site, by="margin")##
## Call:
## adonis(formula = hem.2014 ~ env.2014$tree.height + leaftraits$morphotype + leaftraits$meanSLA + leaftraits$percN + env.2014$z.age + xy.2014$lat + xy.2014$long + env.2014$site, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## env.2014$tree.height 1 0.6339 0.63394 2.37893 0.07631 0.012 *
## leaftraits$morphotype 2 1.0796 0.53981 2.02572 0.12996 0.014 *
## leaftraits$meanSLA 1 0.7449 0.74495 2.79552 0.08968 0.005 **
## leaftraits$percN 1 0.3277 0.32774 1.22988 0.03945 0.251
## env.2014$z.age 1 0.4335 0.43353 1.62688 0.05219 0.099 .
## xy.2014$lat 1 0.7524 0.75241 2.82353 0.09057 0.006 **
## xy.2014$long 1 0.3657 0.36574 1.37249 0.04403 0.208
## env.2014$site 2 0.5050 0.25250 0.94754 0.06079 0.542
## Residuals 13 3.4642 0.26648 0.41702
## Total 23 8.3072 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### try again for all 2014, not just hybrid
env.sub<-env[env$vialID%in%rownames(leaftraits.full),]
hem.sub<-hem[rownames(hem)%in%rownames(leaftraits.full),]
adonis(hem.sub~leaftraits.full$morphotype+leaftraits.full$meanSLA+leaftraits.full$percN+env.sub$z.age+env.sub$site, by="margin")##
## Call:
## adonis(formula = hem.sub ~ leaftraits.full$morphotype + leaftraits.full$meanSLA + leaftraits.full$percN + env.sub$z.age + env.sub$site, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## leaftraits.full$morphotype 2 2.0183 1.00917 3.7114 0.07418 0.001 ***
## leaftraits.full$meanSLA 1 1.2077 1.20774 4.4417 0.04439 0.001 ***
## leaftraits.full$percN 1 0.2891 0.28909 1.0632 0.01062 0.369
## env.sub$z.age 1 0.8705 0.87046 3.2013 0.03199 0.001 ***
## env.sub$site 11 7.5954 0.69049 2.5394 0.27916 0.001 ***
## Residuals 56 15.2271 0.27191 0.55965
## Total 72 27.2081 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1temp<-hem[env$sampleyear=="2014",]
temp<-temp[env$vialID%in%rownames(leaftraits.full),]
mod.leaf<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+leaftraits$meanWatercontent+leaftraits$percP)
mod.leaf.age<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$age)
mod.leaf.tree<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height)
mod.leaf.age.site<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
mod.leaf.site<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
modfull<-cca(hem[rownames(hem)%in%rownames(leaftraits),]~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
m1<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
m2<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+leaftraits$percP+leaftraits$meanWatercontent)
m3<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$tree.height+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
m4<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$z.age+env.2014$site+leaftraits$meanWatercontent+leaftraits$percP)
m5<-cca(hem.2014~leaftraits$percN+leaftraits$morphotype+leaftraits$meanSLA+env.2014$site)
m6<-cca(hem.2014~leaftraits$morphotype+leaftraits$meanSLA+env.2014$site )
anova(m5, by = "margin", step=200 ) ## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = hem.2014 ~ leaftraits$percN + leaftraits$morphotype + leaftraits$meanSLA + env.2014$site)
## Df ChiSquare F Pr(>F)
## leaftraits$percN 1 0.09340 1.0262 0.365
## leaftraits$morphotype 2 0.25951 1.4257 0.044 *
## leaftraits$meanSLA 1 0.11201 1.2307 0.204
## env.2014$site 3 0.48967 1.7934 0.004 **
## Residual 16 1.45618
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 anova(m5,m6)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: hem.2014 ~ leaftraits$percN + leaftraits$morphotype + leaftraits$meanSLA + env.2014$site
## Model 2: hem.2014 ~ leaftraits$morphotype + leaftraits$meanSLA + env.2014$site
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 16 1.4562
## 2 17 1.5496 -1 -0.093399 1.0262 0.36m7<-cca(hem.2014~leaftraits$morphotype+env.2014$site )
anova(m6,m7)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: hem.2014 ~ leaftraits$morphotype + leaftraits$meanSLA + env.2014$site
## Model 2: hem.2014 ~ leaftraits$morphotype + env.2014$site
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 17 1.5496
## 2 18 1.6780 -1 -0.12845 1.4091 0.092 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m8<-cca(hem.2014~env.2014$site )
anova(m7,m8)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: hem.2014 ~ leaftraits$morphotype + env.2014$site
## Model 2: hem.2014 ~ env.2014$site
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 18 1.6780
## 2 20 2.1829 -2 -0.50488 2.7079 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model with just site is bestDo CCA and permanova with data lumped for 2014 and 2015
# subset data in species by site matrix such that it has only 24 rows, one for each site plot combination
# each row should be the sum of collecting in 2014 and 2015
# be careful to exclude anywhere where sampling effort is wildly different
hem$vialID<-rownames(hem)
temp<-melt(hem)## Using vialID as id variables#trial<-dcast(temp, site+plot~variable, fun.aggregate = sum,value.var="value")
trial.env<-merge(temp, env[,1:3], by="vialID")
trial<-dcast(trial.env, site+plot~variable, fun.aggregate = sum,value.var="value")
trial<-trial[trial$site%in%env.both$site,]
hem<-hem[,-84]
rowSums(hem[env$site=="KAIY",])## 25 26 28 30 42 58 141 142 143 163 164 165 166
## 6 38 12 4 4 10 8 22 16 76 18 29 25rowSums(trial[,-c(1:2)])## 1 2 3 4 5 6 7 8 9 10 11 12 25 26 27 28 29 30
## 42 36 130 83 31 12 44 42 65 38 20 18 22 22 12 39 37 136
## 73 74 75 76 77 78
## 55 611 56 592 120 494##################################################
env.trial<-trial[,1:2]
env.trial<-unique(merge(env.trial, env[,c(2,3,10)], by=c("site","plot")))
trial<-trial[,-c(1:2)]
moment<-metaMDS(trial, k=3)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1492291
## Run 1 stress 0.1659777
## Run 2 stress 0.1489775
## ... New best solution
## ... Procrustes: rmse 0.03255228 max resid 0.1093367
## Run 3 stress 0.1661723
## Run 4 stress 0.1489705
## ... New best solution
## ... Procrustes: rmse 0.06705354 max resid 0.1853576
## Run 5 stress 0.148974
## ... Procrustes: rmse 0.0009068024 max resid 0.002178237
## ... Similar to previous best
## Run 6 stress 0.1489704
## ... New best solution
## ... Procrustes: rmse 0.0012843 max resid 0.002781557
## ... Similar to previous best
## Run 7 stress 0.1489776
## ... Procrustes: rmse 0.06670506 max resid 0.1817755
## Run 8 stress 0.1489698
## ... New best solution
## ... Procrustes: rmse 0.001356732 max resid 0.004352775
## ... Similar to previous best
## Run 9 stress 0.1489692
## ... New best solution
## ... Procrustes: rmse 0.001615042 max resid 0.004077342
## ... Similar to previous best
## Run 10 stress 0.1690282
## Run 11 stress 0.1800283
## Run 12 stress 0.1492302
## ... Procrustes: rmse 0.05076003 max resid 0.1856264
## Run 13 stress 0.1489787
## ... Procrustes: rmse 0.06748909 max resid 0.1828223
## Run 14 stress 0.16617
## Run 15 stress 0.1489843
## ... Procrustes: rmse 0.0672218 max resid 0.1824667
## Run 16 stress 0.1489705
## ... Procrustes: rmse 0.0004047275 max resid 0.00093899
## ... Similar to previous best
## Run 17 stress 0.1661743
## Run 18 stress 0.1489707
## ... Procrustes: rmse 0.00210368 max resid 0.004964002
## ... Similar to previous best
## Run 19 stress 0.1492229
## ... Procrustes: rmse 0.05149002 max resid 0.185743
## Run 20 stress 0.1489796
## ... Procrustes: rmse 0.06739322 max resid 0.1828275
## *** Solution reachedmds.fig <- ordiplot(moment, type = "none", main=" all")
text(mds.fig, "species", cex=0.5)
points(mds.fig, "sites", pch = 15, col = "blue", select = env.trial$morphotype == "H")
points(mds.fig, "sites", pch = 19, col = "red", select = env.trial$morphotype == "P")
points(mds.fig, "sites", pch = 19, col = "green", select = env.trial$morphotype == "G")
adonis(trial~env.trial$morphotype, by="margin")##
## Call:
## adonis(formula = trial ~ env.trial$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## env.trial$morphotype 2 1.3344 0.66719 2.4863 0.19146 0.001 ***
## Residuals 21 5.6352 0.26834 0.80854
## Total 23 6.9696 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# with 2014 and 2015 data summed pretty similar results to when with twice as many data points
# pretend it has vialID
env.trial<-merge(env.trial, env.2014[,1:3], by=c("site","plot"))
rownames(trial)<-env.trial$vialID
# make sure all same order
env.trial<-env.trial[order(env.trial$vialID),]
trial<-trial[order(rownames(trial)),]
m1<-adonis(trial~leaftraits$meanWatercontent+env.trial$morphotype+leaftraits$meanSLA+leaftraits$percN+leaftraits$percP, by="margin")
# even as first variable water content is not sig
adonis(trial~env.trial$morphotype+leaftraits$meanSLA+leaftraits$percN+leaftraits$percP, by="margin")##
## Call:
## adonis(formula = trial ~ env.trial$morphotype + leaftraits$meanSLA + leaftraits$percN + leaftraits$percP, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## env.trial$morphotype 2 1.3369 0.66845 2.73289 0.19182 0.001 ***
## leaftraits$meanSLA 1 0.5869 0.58685 2.39928 0.08420 0.014 *
## leaftraits$percN 1 0.2387 0.23869 0.97584 0.03425 0.470
## leaftraits$percP 1 0.4044 0.40440 1.65336 0.05802 0.085 .
## Residuals 18 4.4027 0.24460 0.63171
## Total 23 6.9696 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# percN NS in third place, what about up front
adonis(trial~leaftraits$percN+leaftraits$meanSLA+env.trial$morphotype+leaftraits$percP, by="margin")##
## Call:
## adonis(formula = trial ~ leaftraits$percN + leaftraits$meanSLA + env.trial$morphotype + leaftraits$percP, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## leaftraits$percN 1 0.7208 0.72078 2.9468 0.10342 0.002 **
## leaftraits$meanSLA 1 0.3543 0.35427 1.4484 0.05083 0.156
## env.trial$morphotype 2 1.0874 0.54370 2.2229 0.15602 0.004 **
## leaftraits$percP 1 0.4044 0.40440 1.6534 0.05802 0.092 .
## Residuals 18 4.4027 0.24460 0.63171
## Total 23 6.9696 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# does better in first place, seems to correlate alot with SLA. what if only including SLA- sig in last place?
adonis(trial~env.trial$morphotype+leaftraits$percP+leaftraits$meanSLA, by="margin")##
## Call:
## adonis(formula = trial ~ env.trial$morphotype + leaftraits$percP + leaftraits$meanSLA, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## env.trial$morphotype 2 1.3369 0.66845 2.6954 0.19182 0.002 **
## leaftraits$percP 1 0.2545 0.25451 1.0263 0.03652 0.425
## leaftraits$meanSLA 1 0.6662 0.66618 2.6862 0.09558 0.004 **
## Residuals 19 4.7120 0.24800 0.67608
## Total 23 6.9696 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# SLA still sig. percP is not in second place (it is sig in first place)
adonis(trial~leaftraits$meanSLA+env.trial$morphotype, by="margin")##
## Call:
## adonis(formula = trial ~ leaftraits$meanSLA + env.trial$morphotype, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## leaftraits$meanSLA 1 0.8404 0.84042 3.3311 0.12058 0.002 **
## env.trial$morphotype 2 1.0833 0.54167 2.1470 0.15544 0.013 *
## Residuals 20 5.0458 0.25229 0.72398
## Total 23 6.9696 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# mean SLA and morphotype sig
# regardless of which comes first, together R2 0.28, in both cases ~p<0.01
# what about xy:
adonis(trial~leaftraits$percN+leaftraits$percP+env.trial$morphotype+leaftraits$meanSLA+xy.2014$lat, by="margin")##
## Call:
## adonis(formula = trial ~ leaftraits$percN + leaftraits$percP + env.trial$morphotype + leaftraits$meanSLA + xy.2014$lat, by = "margin")
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## leaftraits$percN 1 0.7208 0.72078 3.1357 0.10342 0.002 **
## leaftraits$percP 1 0.5768 0.57683 2.5095 0.08276 0.010 **
## env.trial$morphotype 2 0.8222 0.41112 1.7886 0.11798 0.037 *
## leaftraits$meanSLA 1 0.4470 0.44700 1.9446 0.06414 0.024 *
## xy.2014$lat 1 0.4951 0.49509 2.1539 0.07104 0.019 *
## Residuals 17 3.9076 0.22986 0.56067
## Total 23 6.9696 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod4<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$lat+xy.2014$long)
mod3<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent)
mod5<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$lat)
anova(mod3,mod4,mod5)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 3: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 18 1.9822
## 2 15 1.5024 3 0.47982 1.5968 0.015 *
## 3 16 1.6095 -1 -0.10710 1.0693 0.377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(mod3,mod4)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 18 1.9822
## 2 15 1.5024 3 0.47982 1.5968 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model shortest is best
# try removing latitude
mod6<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$long)
anova(mod4,mod6)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$long
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 15 1.5024
## 2 16 1.5896 -1 -0.087169 0.8703 0.642# shortest model is best
mod3<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent)
mod2<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN)
anova(mod2,mod3)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 19 2.2269
## 2 18 1.9822 1 0.24468 2.2219 0.007 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model 3 seems better
mod<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$meanWatercontent)
anova(mod,mod3)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 19 2.1035
## 2 18 1.9822 1 0.12131 1.1016 0.357# no difference
mod6<-cca(decostand(trial, "total")~leaftraits$meanSLA+leaftraits$meanWatercontent)
anova(mod,mod6)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + leaftraits$meanWatercontent
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 19 2.1035
## 2 21 2.5572 -2 -0.45369 2.049 0.002 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod7<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype)
anova(mod,mod7)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$meanWatercontent
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 19 2.1035
## 2 20 2.3389 -1 -0.23535 2.1258 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod8<-cca(decostand(trial, "total")~leaftraits$meanSLA)
anova(mod4,mod7,mod8)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype
## Model 3: decostand(trial, "total") ~ leaftraits$meanSLA
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 15 1.5024
## 2 20 2.3389 -5 -0.83648 1.6703 0.002 **
## 3 22 2.7205 -2 -0.38163 1.9051 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod9<-cca(decostand(trial, "total")~env.trial$morphotype)
anova(mod4,mod7,mod9)## Permutation tests for cca under reduced model
## Permutation: free
## Number of permutations: 999
##
## Model 1: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long
## Model 2: decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype
## Model 3: decostand(trial, "total") ~ env.trial$morphotype
## ResDf ResChiSquare Df ChiSquare F Pr(>F)
## 1 15 1.5024
## 2 20 2.3389 -5 -0.83648 1.6703 0.001 ***
## 3 21 2.5950 -1 -0.25611 2.5571 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# can't compare model with just morphotype or just SLA, but F value in the same comparison (full model, model with the two variables and model with just one) is higher for morphotype
# best model mod9
mod4<-cca(decostand(trial, "total")~leaftraits$meanSLA+env.trial$morphotype+leaftraits$percN+leaftraits$meanWatercontent+leaftraits$percP+xy.2014$lat+xy.2014$long)
plot(mod4, display = "sites", scaling = 3, type="points")
text(mod4, scaling = 3, display = "bp", pch=4)
text(mod4, display = "species", scaling = 3, col="red", cex=0.5)
anova(mod4, by = "margin", step=200 )## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = decostand(trial, "total") ~ leaftraits$meanSLA + env.trial$morphotype + leaftraits$percN + leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat + xy.2014$long)
## Df ChiSquare F Pr(>F)
## leaftraits$meanSLA 1 0.14319 1.4297 0.146
## env.trial$morphotype 2 0.30895 1.5423 0.037 *
## leaftraits$percN 1 0.11347 1.1329 0.299
## leaftraits$meanWatercontent 1 0.14256 1.4233 0.124
## leaftraits$percP 1 0.11941 1.1922 0.303
## xy.2014$lat 1 0.08717 0.8703 0.633
## xy.2014$long 1 0.10710 1.0693 0.397
## Residual 15 1.50239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod4## Call: cca(formula = decostand(trial, "total") ~ leaftraits$meanSLA
## + env.trial$morphotype + leaftraits$percN +
## leaftraits$meanWatercontent + leaftraits$percP + xy.2014$lat +
## xy.2014$long)
##
## Inertia Proportion Rank
## Total 3.0213 1.0000
## Constrained 1.5189 0.5027 8
## Unconstrained 1.5024 0.4973 15
## Inertia is scaled Chi-square
## 38 species (variables) deleted due to missingness
##
## Eigenvalues for constrained axes:
## CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 CCA8
## 0.4530 0.3274 0.2532 0.2008 0.1061 0.0760 0.0602 0.0422
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8 CA9
## 0.31373 0.27569 0.16823 0.16084 0.13136 0.09011 0.08342 0.06841 0.05593
## CA10 CA11 CA12 CA13 CA14 CA15
## 0.04687 0.03611 0.02644 0.01971 0.01593 0.00962# morphotype, mean water and percP are sig, SLA marginally
# on figure, percN, percP and SLA completely overlaying eachother. morphotype H and P also fall in same spot
# explains 41% variation
# explore:
# seems like all of the dimensions around leaftraits are captured by meanSLA (cor strongly with percN, percP), morphotype and water content
vare.cca <- cca(decostand(trial, "total") ~ env.trial$morphotype)
plot(vare.cca, display = "sites", scaling = 3, type="points")
text(vare.cca, scaling = 3, display = "bp", pch=4)
text(vare.cca, display = "species", scaling = 3, col="red", cex=0.5)
anova(vare.cca, by = "margin", step=200 )## Permutation test for cca under reduced model
## Marginal effects of terms
## Permutation: free
## Number of permutations: 999
##
## Model: cca(formula = decostand(trial, "total") ~ env.trial$morphotype)
## Df ChiSquare F Pr(>F)
## env.trial$morphotype 2 0.42632 1.725 0.001 ***
## Residual 21 2.59498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1vare.cca## Call: cca(formula = decostand(trial, "total") ~
## env.trial$morphotype)
##
## Inertia Proportion Rank
## Total 3.0213 1.0000
## Constrained 0.4263 0.1411 2
## Unconstrained 2.5950 0.8589 21
## Inertia is scaled Chi-square
## 38 species (variables) deleted due to missingness
##
## Eigenvalues for constrained axes:
## CCA1 CCA2
## 0.23183 0.19449
##
## Eigenvalues for unconstrained axes:
## CA1 CA2 CA3 CA4 CA5 CA6 CA7 CA8
## 0.5219 0.3710 0.2926 0.2466 0.2218 0.1693 0.1407 0.1262
## (Showed only 8 of all 21 unconstrained eigenvalues)# together 30% variation(iii) Which insect traits account for herbivore community responses to host-plant phenotype
Are communities distinctly different between different morphotype trees across all youngish sites?
moment<-metaMDS(hem.both, k=3)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1790701
## Run 1 stress 0.1860155
## Run 2 stress 0.1839765
## Run 3 stress 0.1892048
## Run 4 stress 0.1865272
## Run 5 stress 0.1791475
## ... Procrustes: rmse 0.008537058 max resid 0.04929761
## Run 6 stress 0.1832359
## Run 7 stress 0.1791218
## ... Procrustes: rmse 0.01407452 max resid 0.08428514
## Run 8 stress 0.1889459
## Run 9 stress 0.179072
## ... Procrustes: rmse 0.0007812932 max resid 0.002826476
## ... Similar to previous best
## Run 10 stress 0.1790731
## ... Procrustes: rmse 0.0008433542 max resid 0.00345565
## ... Similar to previous best
## Run 11 stress 0.1832108
## Run 12 stress 0.1826911
## Run 13 stress 0.1860141
## Run 14 stress 0.1854629
## Run 15 stress 0.186774
## Run 16 stress 0.1790705
## ... Procrustes: rmse 0.0008900677 max resid 0.004086788
## ... Similar to previous best
## Run 17 stress 0.1791349
## ... Procrustes: rmse 0.01539417 max resid 0.08989465
## Run 18 stress 0.1805679
## Run 19 stress 0.1791347
## ... Procrustes: rmse 0.01156497 max resid 0.06346796
## Run 20 stress 0.1791222
## ... Procrustes: rmse 0.01373238 max resid 0.08235413
## *** Solution reachedmds.fig <- ordiplot(moment, type = "none", main="all hybrid zone sites")
#text(mds.fig, "species")
points(mds.fig, "sites", pch = 19, col = "green", select = env.both$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env.both$morphotype == "H")
points(mds.fig, "sites", pch = 17, col = "red", select = env.both$morphotype == "P")
Lab<-c("Glabrous", "Hybrid", "Pubescent")
color<-c("green", "blue", "red")
legend( "topright", legend=Lab,pch=c(19,15,17),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)
## across all sites (n=170)
moment<-metaMDS(hem, k=3)## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2143791
## Run 1 stress 0.2144895
## ... Procrustes: rmse 0.02158443 max resid 0.09391575
## Run 2 stress 0.2143115
## ... New best solution
## ... Procrustes: rmse 0.04348152 max resid 0.1854511
## Run 3 stress 0.2143784
## ... Procrustes: rmse 0.03803577 max resid 0.1705378
## Run 4 stress 0.2147627
## ... Procrustes: rmse 0.02768923 max resid 0.1083661
## Run 5 stress 0.2148255
## Run 6 stress 0.2135422
## ... New best solution
## ... Procrustes: rmse 0.03912804 max resid 0.1809095
## Run 7 stress 0.2141577
## Run 8 stress 0.2132945
## ... New best solution
## ... Procrustes: rmse 0.01515823 max resid 0.08904171
## Run 9 stress 0.2163937
## Run 10 stress 0.2144438
## Run 11 stress 0.2133151
## ... Procrustes: rmse 0.002346661 max resid 0.0125573
## Run 12 stress 0.2142829
## Run 13 stress 0.2134329
## ... Procrustes: rmse 0.01233698 max resid 0.1037661
## Run 14 stress 0.2135966
## ... Procrustes: rmse 0.01799393 max resid 0.08714684
## Run 15 stress 0.2148402
## Run 16 stress 0.2134204
## ... Procrustes: rmse 0.01341335 max resid 0.07907385
## Run 17 stress 0.2146835
## Run 18 stress 0.2144819
## Run 19 stress 0.2180064
## Run 20 stress 0.2142619
## *** No convergence -- monoMDS stopping criteria:
## 6: no. of iterations >= maxit
## 14: stress ratio > sratmaxmds.fig <- ordiplot(moment, type = "none", main="all sites")
text(mds.fig, "species",cex=0.5)
points(mds.fig, "sites", pch = 19, col = "grey", select = env$morphotype == "G")
points(mds.fig, "sites", pch = 15, col = "blue", select = env$morphotype == "H")
points(mds.fig, "sites", pch = 17, col = "red", select = env$morphotype == "P")
Lab<-c("Glabrous", "Hybrid", "Pubescent")
color<-c("grey", "blue", "red")
legend( "topright", legend=Lab,pch=c(19,15,17),pt.bg=c(color),pt.cex=1.1,cex=0.9,xpd=TRUE)
# for all sites without double measures (n=97)
#moment<-metaMDS(hem.all, k=3)
#mds.fig <- ordiplot(moment, type = "none", main="all sites")
#text(mds.fig, "species",cex=0.5)
#points(mds.fig, "sites", pch = 19, col = "grey", select = env$morphotype == "G")
#points(mds.fig, "sites", pch = 15, col = "blue", select = env$morphotype == "H")
#points(mds.fig, "sites", pch = 17, col = "red", select = env$morphotype == "P")phosphorus and morphotype are only sig when only/first term
mean SLA and mean water content are always sig no matter their placement
perc N significant unless after SLA and before water content
What does diversity and abundance look like across these sites with multiple phenotypes?
plot(1/(1-ChaoSimpson(as.data.frame(t(hem.ero)) )[,2])~env[env$site=="ERO",]$morphotype)
plot(1/(1-ChaoSimpson(as.data.frame(t(hem.kaiy)) )[,2])~env[env$site=="KAIY",]$morphotype)## Warning in bplt(at[i], wid = width[i], stats = z$stats[, i], out =
## z$out[z$group == : Outlier (Inf) in boxplot 4 is not drawn
plot(1/(1-ChaoSimpson(as.data.frame(t(hem.olaa)) )[,2])~env[env$site=="OLAA",]$morphotype)
plot(1/(1-ChaoSimpson(as.data.frame(t(hem.ali)) )[,2])~env[env$site=="ALI",]$morphotype)## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
plot(ChaoRichness(as.data.frame(t(hem.ero)) )[,2]~env[env$site=="ERO",]$morphotype)
plot(ChaoRichness(as.data.frame(t(hem.kaiy)) )[,2]~env[env$site=="KAIY",]$morphotype)
plot(ChaoRichness(as.data.frame(t(hem.olaa)) )[,2]~env[env$site=="OLAA",]$morphotype)
plot(ChaoRichness(as.data.frame(t(hem.ali)) )[,2]~env[env$site=="ALI",]$morphotype)
plot(rowSums(hem.ero)~env[env$site=="ERO",]$morphotype)
plot(rowSums(hem.kaiy)~env[env$site=="KAIY",]$morphotype)
plot(rowSums(hem.olaa)~env[env$site=="OLAA",]$morphotype)
plot(rowSums(hem.ali)~env[env$site=="ALI",]$morphotype)
anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.ero)) )[,2])~env[env$site=="ERO",]$morphotype))## Analysis of Variance Table
##
## Response: 1/(1 - ChaoSimpson(as.data.frame(t(hem.ero)))[, 2])
## Df Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "ERO", ]$morphotype 1 9.503 9.5031 0.9171 0.3608
## Residuals 10 103.626 10.3626#anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.kaiy)) )[,2])~env[env$site=="KAIY",]$morphotype))
anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.olaa)) )[,2])~env[env$site=="OLAA",]$morphotype))## Analysis of Variance Table
##
## Response: 1/(1 - ChaoSimpson(as.data.frame(t(hem.olaa)))[, 2])
## Df Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "OLAA", ]$morphotype 1 0.2931 0.2931 0.1316 0.7244
## Residuals 10 22.2796 2.2280anova(lm(1/(1-ChaoSimpson(as.data.frame(t(hem.ali)) )[,2])~env[env$site=="ALI",]$morphotype))## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.## Analysis of Variance Table
##
## Response: 1/(1 - ChaoSimpson(as.data.frame(t(hem.ali)))[, 2])
## Df Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "ALI", ]$morphotype 1 46.116 46.116 1.8435 0.2044
## Residuals 10 250.158 25.016anova(lm(rowSums(hem.kaiy)~env[env$site=="KAIY",]$morphotype))## Analysis of Variance Table
##
## Response: rowSums(hem.kaiy)
## Df Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "KAIY", ]$morphotype 1 1589.0 1589.00 5.803 0.03468 *
## Residuals 11 3012.1 273.83
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# p<0.05
anova(lm(rowSums(hem.olaa)~env[env$site=="OLAA",]$morphotype))## Analysis of Variance Table
##
## Response: rowSums(hem.olaa)
## Df Sum Sq Mean Sq F value Pr(>F)
## env[env$site == "OLAA", ]$morphotype 1 81667 81667 4.6195 0.05713 .
## Residuals 10 176786 17679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# p<0.1
# Across all 12 sites:
env.both$morphotype<-factor(env.both$morphotype, levels=c("H", "G", "P"))
plot(rowSums(hem.both)~env.both$morphotype, ylim=c(0,300), ylab="Abundance")
anova(lm(rowSums(hem.both)~env.both$morphotype))## Analysis of Variance Table
##
## Response: rowSums(hem.both)
## Df Sum Sq Mean Sq F value Pr(>F)
## env.both$morphotype 2 91547 45773 5.9593 0.004996 **
## Residuals 46 353328 7681
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(aov(rowSums(hem.both)~env.both$morphotype))## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = rowSums(hem.both) ~ env.both$morphotype)
##
## $`env.both$morphotype`
## diff lwr upr p adj
## G-H 82.796154 7.178384 158.41392 0.0289491
## P-H -8.841346 -88.095293 70.41260 0.9605958
## P-G -91.637500 -162.829276 -20.44572 0.0086619# p<0.05
tem<-1-ChaoSimpson(as.data.frame(t(hem.both), na.rm=TRUE) )## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.#tem[10,2]<-1
#tem[65,2]<-1
tem<-ifelse(tem[,2]==0,0,1/tem[,2])
anova(lm(tem~env.both$morphotype))## Analysis of Variance Table
##
## Response: tem
## Df Sum Sq Mean Sq F value Pr(>F)
## env.both$morphotype 2 68.32 34.160 2.6805 0.07922 .
## Residuals 46 586.22 12.744
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lmer(tem~env.both$morphotype+(1|env.both$site)))## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## env.both$morphotype 2 53.867 26.933 2.2435tem<-exp(ChaoShannon(as.data.frame(t(hem.both)) )[,2])## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.#tem[2]<-1
#tem[26]<-1
anova(lm(tem~env.both$morphotype))## Analysis of Variance Table
##
## Response: tem
## Df Sum Sq Mean Sq F value Pr(>F)
## env.both$morphotype 2 68.43 34.217 3.074 0.05584 .
## Residuals 46 512.04 11.131
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lmer(tem~env.both$morphotype+(1|env.both$site)))## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## env.both$morphotype 2 66.67 33.335 3.0035tem<-ChaoRichness(as.data.frame(t(hem.both)))[,2]
anova(lm(tem~env.both$morphotype))## Analysis of Variance Table
##
## Response: tem
## Df Sum Sq Mean Sq F value Pr(>F)
## env.both$morphotype 2 233.82 116.912 2.6616 0.08058 .
## Residuals 46 2020.60 43.926
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(lmer(tem~env.both$morphotype+(1|env.both$site)))## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## env.both$morphotype 2 241.78 120.89 2.7891plot(tem~env.both$morphotype, ylab="Species Richness")
# richness does diff sig between morphotype, diversity shannon or simpson does not.
# ie the difference between glabrous and hybrid and glabrous and pubescent is in rare species- glabrous has more rare species than the other two morphotypes, but once you weigh rare species less, there are no differences in diversity
### plots for all sites
plot(diversity(hem)/log(specnumber(hem))~env$morphotype, ylab="Pielou's Evenness")
plot(rowSums(hem)~env$morphotype, ylab="Abundance",xlab="Tree morphotype", main="across big island (n=170)", ylim=c(0,350))
tem<-exp(ChaoShannon(as.data.frame(t(hem)) )[,2])## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.
## Warning in BootstrapFun.abun(x = x, FunName, datatype, B): This site has
## only one species. Estimation is not robust.plot(tem~env$morphotype, ylab="Exp Shannon Index",main="across big island (n=170)", xlab="Tree morphotype")
How much variation is explained by spatial proximity?
GLMER
# dataset with 2014 and 2015 data
hem.both$vialID<-rownames(hem.both)
hem.both.long<-melt(hem.both)## Using vialID as id variableshem.both.long<-merge(env.both, hem.both.long, by="vialID")
leaftraits$vialID<-rownames(leaftraits)
temp<-merge(leaftraits,env[,c(1:3)], by="vialID")
temp<-temp[,-1]
hem.both.long<-hem.both.long[,-1]
hem.both.long<-merge(hem.both.long, temp, by=c("site", "plot"))
hem.both.long$percN.standard<-(hem.both.long$percN-mean(hem.both.long$percN))/sd(hem.both.long$percN)
hem.both.long$percP.standard<-(hem.both.long$percP-mean(hem.both.long$percP))/sd(hem.both.long$percP)
hem.both.long$meanSLA.standard<-(hem.both.long$meanSLA-mean(hem.both.long$meanSLA))/sd(hem.both.long$meanSLA)
hem.both.long$meanWatercontent.standard<-(hem.both.long$meanWatercontent-mean(hem.both.long$meanWatercontent))/sd(hem.both.long$meanWatercontent)
############## analysis data 2014, 2015
modelP=glmer(value ~morphotype.x +meanSLA.standard +percP.standard+(1+percP.standard|variable)+(1|site), data =hem.both.long, family =poisson)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00111975
## (tol = 0.001, component 1)qqnorm(resid(modelP))
# overdispersed? negative binomial more appropriate
# best model is modelP, including morphotype, SLA and P, with species specific response to P and site also as random effect
modelSLA=glmer.nb(value ~meanSLA.standard+percP.standard+(1+meanSLA.standard|variable)+(1|site), data =hem.both.long)
modelP=glmer.nb(value ~meanSLA.standard+percP.standard+(1+percP.standard|variable)+(1|site), data =hem.both.long)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0279573
## (tol = 0.001, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0294294
## (tol = 0.001, component 1)modelN=glmer.nb(value ~percN.standard+percP.standard+(1+percN.standard|variable)+(1|site), data =hem.both.long)
anova(modelP,modelSLA,modelN)## Data: hem.both.long
## Models:
## modelP: value ~ meanSLA.standard + percP.standard + (1 + percP.standard |
## modelP: variable) + (1 | site)
## modelSLA: value ~ meanSLA.standard + percP.standard + (1 + meanSLA.standard |
## modelSLA: variable) + (1 | site)
## modelN: value ~ percN.standard + percP.standard + (1 + percN.standard |
## modelN: variable) + (1 | site)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modelP 8 3009.1 3054.5 -1496.6 2993.1
## modelSLA 8 3011.2 3056.6 -1497.6 2995.2 0 0 1
## modelN 8 3015.2 3060.5 -1499.6 2999.2 0 0 1modelPalone=glmer.nb(value ~percP.standard+(1+percP.standard|variable)+(1|site), data =hem.both.long)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0289687
## (tol = 0.001, component 1)modelPmorph=glmer.nb(value ~morphotype.x+percP.standard+(1+percP.standard|variable)+(1|site), data =hem.both.long)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0193233
## (tol = 0.001, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge: degenerate Hessian with 1
## negative eigenvaluesmodelPSLA=glmer.nb(value ~meanSLA.standard+(1+percP.standard|variable)+(1|site), data =hem.both.long)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge: degenerate Hessian with 1
## negative eigenvaluesanova(modelP,modelPalone,modelPmorph)## Data: hem.both.long
## Models:
## modelPalone: value ~ percP.standard + (1 + percP.standard | variable) + (1 |
## modelPalone: site)
## modelP: value ~ meanSLA.standard + percP.standard + (1 + percP.standard |
## modelP: variable) + (1 | site)
## modelPmorph: value ~ morphotype.x + percP.standard + (1 + percP.standard |
## modelPmorph: variable) + (1 | site)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modelPalone 7 3011.3 3051.0 -1498.7 2997.3
## modelP 8 3009.1 3054.5 -1496.6 2993.1 4.1927 1 0.0405987
## modelPmorph 9 2999.5 3050.5 -1490.7 2981.5 11.6537 1 0.0006407
##
## modelPalone
## modelP *
## modelPmorph ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1modelwP=glmer.nb(value ~morphotype.x+percP.standard+(1|variable)+(1|site), data =hem.both.long)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0253589
## (tol = 0.001, component 1)modelPwallP=glmer.nb(value ~morphotype.x+(1|variable)+(1|site), data =hem.both.long)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0224301
## (tol = 0.001, component 1)anova(modelPmorph,modelwP,modelPwallP)## Data: hem.both.long
## Models:
## modelPwallP: value ~ morphotype.x + (1 | variable) + (1 | site)
## modelwP: value ~ morphotype.x + percP.standard + (1 | variable) + (1 |
## modelwP: site)
## modelPmorph: value ~ morphotype.x + percP.standard + (1 + percP.standard |
## modelPmorph: variable) + (1 | site)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## modelPwallP 6 3009.2 3043.2 -1498.6 2997.2
## modelwP 7 3011.1 3050.8 -1498.6 2997.1 0.0485 1 0.8257047
## modelPmorph 9 2999.5 3050.5 -1490.7 2981.5 15.6670 2 0.0003962
##
## modelPwallP
## modelwP
## modelPmorph ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# model with morphotype and perc P, and random effect of percP is best model
summary(modelPmorph)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.2625) ( log )
## Formula: value ~ morphotype.x + percP.standard + (1 + percP.standard |
## variable) + (1 | site)
## Data: hem.both.long
##
## AIC BIC logLik deviance df.resid
## 2999.5 3050.5 -1490.7 2981.5 2133
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5098 -0.3360 -0.2173 -0.1338 11.8310
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## variable (Intercept) 3.5283 1.8784
## percP.standard 0.2258 0.4752 0.14
## site (Intercept) 0.1920 0.4382
## Number of obs: 2142, groups: variable, 42; site, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.51629 0.40967 -6.142 8.13e-10 ***
## morphotype.xG 0.88645 0.22141 4.004 6.24e-05 ***
## morphotype.xP 0.10121 0.28266 0.358 0.72
## percP.standard -0.01331 0.17780 -0.075 0.94
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mrph.G mrph.P
## morphtyp.xG -0.305
## morphtyp.xP -0.235 0.126
## prcP.stndrd 0.043 -0.394 0.337dotplot(ranef(modelP, condVar=T))## $variable
##
## $site
# check for overdispersion:
pcsq<-sum(resid(modelP,type="pearson")^2)
rdf<-df.residual(modelP)
pchisq(pcsq,rdf,lower.tail=F)## [1] 1pcsq/rdf## [1] 0.8416044library(blmeco)## Loading required package: MASSdispersion_glmer(modelP)## [1] 0.6364843# alternatively could use package DHArma
################# Figures 2015 and 2014 data
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000", "green")
ggplot(data=hem.both.long, aes(x=percP.standard, y=log(value), group=variable, colour=variable)) +
# geom_point()+
geom_smooth(method=lm, se=F,level=0.95)+
ylab("Abundance")+ xlab("standardized perc P")+
theme_bw()## Warning: Removed 1795 rows containing non-finite values (stat_smooth).
hemtop10<-hem.both.long[hem.both.long$variable==c("GRE_PSI"),]
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("OCE_VUL"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("OPU_SP"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("LEI_SP"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("ORT_MET"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("KOA_HAW"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("SAR_ADO"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("PAR_PYR"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("TYM_TYM"),])
hemtop10<-rbind(hemtop10,hem.both.long[hem.both.long$variable==c("OLI_SP"),])
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000", "green", "#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7", "#000000", "green")
formatlines<-c(1,2,3,4,5,6,7,1,2,3,4,5,6,7,1,2,3,4,5,6)
ggplot(data=hemtop10, aes(x=percP.standard, y=log(value), group=variable, linetype=variable, colour=variable)) +
geom_smooth(method=lm, se=F,level=0.95)+
ylab("Abundance")+ xlab("standardized foliar phosphorus %")+
scale_linetype_manual(values=formatlines, name="Species")+
scale_colour_manual(values=cbPalette, name="Species")+
theme_bw()## Warning: Removed 297 rows containing non-finite values (stat_smooth).
ggplot(data=hemtop10, aes(x=meanSLA.standard, y=log(value), group=variable, linetype=variable, colour=variable)) +
geom_smooth(method=lm, se=F,level=0.95)+
ylab("Abundance")+ xlab("standardized SLA%")+
scale_linetype_manual(values=formatlines, name="Species")+
scale_colour_manual(values=cbPalette, name="Species")+
theme_bw()## Warning: Removed 297 rows containing non-finite values (stat_smooth).
ggplot(data=hemtop10, aes(x=morphotype.x, y=log(value), group=variable, linetype=variable, colour=variable)) +
geom_point()+
ylab("Abundance")+ xlab("standardized foliar phosphorus %")+
scale_shape_manual(values=formatlines, name="Species")+
scale_colour_manual(values=cbPalette, name="Species")+
theme_bw()
Conclusions
(i) How do herbivore communities respond to phenotypic variation?
There is variation in species composition across sites and metrosideros traits. NMDS plots indicate differences between sites and volcanoes in species composition.
(ii) Which plant traits account for herbivore community responses to host-plant phenotype?
CCA results point to the importance of foliar nitrogen, leaf morphotype and specific leaf area. However, there is still a large part of variation not explained.
(iii) Which insect traits account for herbivore community responses to host-plant phenotype
Community composition varies across metrosideros traits. Feeding guild is correlated with foliar %N, where leaf chewers (sig) and sap feeders (NS) are more abundant on higher nutrient trees. Feeding guild is correlated with morphotype, where high trichome density morphotype is correlated with higher abundance gallers and seed feeders, and lower abundances of chewers and sap feeders. Morphotype is positively correlated with proportion specialists (pubescent leaves show relatively higher abundances of specialists). Location of the nymph varies little across morphotypes,